Based on basic code for calculating and mapping climatic novelty using the sigma dissimilarity metric. * written by Colin Mahony * Supp Material of “A closer look at novel climates: new methods and insights at continental to landscape scales”
This code uses the annnual climate normals and monthly ICV for the novelty calculation
##################################
#### System setup ####
##################################
## Specify location of data ####
setwd("/Users/lotterhos/Documents/GitHub/OceanClimateNovelty")
#setwd("~/Desktop/PostDoc/OceanClimateNovelty/")
source("src/0-Novelty_Oceans_Functions.R")
## Warning: package 'raster' was built under R version 3.6.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.6.2
## Loading required package: ade4
## Loading required package: adehabitatMA
## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp
##
## Attaching package: 'adehabitatMA'
## The following object is masked from 'package:raster':
##
## buffer
## Loading required package: CircStats
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following objects are masked from 'package:raster':
##
## area, select
## Loading required package: boot
##
## Attaching package: 'data.table'
## The following object is masked from 'package:raster':
##
## shift
## ── Attaching packages ────────────────────── tidyverse 1.2.1 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0
## Warning: package 'ggplot2' was built under R version 3.6.2
## ── Conflicts ───────────────────────── tidyverse_conflicts() ──
## x dplyr::between() masks data.table::between()
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks data.table::first()
## x dplyr::id() masks adehabitatLT::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks data.table::last()
## x dplyr::select() masks MASS::select(), raster::select()
## x purrr::transpose() masks data.table::transpose()
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.4-0 (2019-11-01) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: maps
##
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
##
## map
## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code
## rgdal: version: 1.4-7, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1
## Warning: replacing previous import 'sf::st_make_valid' by
## 'lwgeom::st_make_valid' when loading 'tmap'
## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo
## Checking rgeos availability: TRUE
## Warning: package 'sf' was built under R version 3.6.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
##
## Attaching package: 'fasterize'
## The following object is masked from 'package:graphics':
##
## plot
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
## The following objects are masked from 'package:data.table':
##
## dcast, melt
##################################
#### Read in the input data ####
##################################
path <- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Data/Lotterhos/"
path2<- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Data/"
dat1 <- fread(paste0(path,"Katie_Temp_Arag_1800_2000.txt"), sep = ",")
dat4.5 <- fread(paste0(path,"Katie_Temp_Arag_2070_2100_RCP45.txt"), sep = ",")
dat8.5 <- fread(paste0(path,"Katie_Temp_Arag_2070_2100_RCP85.txt"), sep = ",")
ICV <- fread(paste0(path2, "ICV_Temp_Arag_pH_1960_2020.txt"))
results_path <- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Results/"
results_dir <- paste0(results_path, "HemisphereRestricted/")
head(dat1)
## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 1800 1 -0.05589 1.8167 8.0967
## 2: 1 20.5 -69.5 1800 2 -0.49722 2.0572 8.1706
## 3: 1 20.5 -69.5 1800 3 -0.88001 1.9321 8.1444
## 4: 1 20.5 -69.5 1800 4 -1.24800 1.8865 8.1396
## 5: 1 20.5 -69.5 1800 5 -2.23610 1.8170 8.1375
## 6: 1 20.5 -69.5 1800 6 -0.21900 1.9615 8.1438
head(dat8.5)
## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 2070 1 0.503480 0.95276 7.7819
## 2: 1 20.5 -69.5 2070 2 0.062159 1.00710 7.8161
## 3: 1 20.5 -69.5 2070 3 -0.320630 0.96601 7.8035
## 4: 1 20.5 -69.5 2070 4 -0.688590 0.94601 7.8008
## 5: 1 20.5 -69.5 2070 5 -1.676700 0.90633 7.7981
## 6: 1 20.5 -69.5 2070 6 0.340380 0.98976 7.8048
head(dat4.5)
## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 2070 1 0.58164 1.1518 7.8685
## 2: 1 20.5 -69.5 2070 2 0.14032 1.2362 7.9107
## 3: 1 20.5 -69.5 2070 3 -0.24248 1.1804 7.8955
## 4: 1 20.5 -69.5 2070 4 -0.61043 1.1555 7.8924
## 5: 1 20.5 -69.5 2070 5 -1.59860 1.1084 7.8899
## 6: 1 20.5 -69.5 2070 6 0.41854 1.2072 7.8965
head(ICV)
## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 1960 1 -0.20420 1.5702 8.0291
## 2: 1 20.5 -69.5 1960 2 -0.64552 1.7479 8.0921
## 3: 1 20.5 -69.5 1960 3 -1.02830 1.6505 8.0698
## 4: 1 20.5 -69.5 1960 4 -1.39630 1.6131 8.0657
## 5: 1 20.5 -69.5 1960 5 -2.38440 1.5515 8.0635
## 6: 1 20.5 -69.5 1960 6 -0.36730 1.6803 8.0699
unique(dat1$Year)
## [1] 1800 1810 1820 1830 1970 1980 1990 2000
unique(dat8.5$Year)
## [1] 2070 2080 2090 2100
unique(dat4.5$Year)
## [1] 2070 2080 2090 2100
unique(ICV$Year)
## [1] 1960 1970 1980 1990 2000 2010 2020
hist(c(dat1$Arag, dat8.5$Arag))
hist(log10(c(dat1$Arag, dat8.5$Arag)))
pdf(paste0(results_dir, "ICV_Lat.pdf"), width=9, height=7)
par(mfrow=c(3,1), mar=c(4,5,0.5,0.5), oma=c(2,0,0,0), las=0)
lat <- sort(unique(ICV$Lat))
length(lat)
## [1] 169
lat[1]
## [1] -78.5
lat[42]
## [1] -37.5
lat[80]
## [1] 0.5
lat[117]
## [1] 37.5
lat[158]
## [1] 78.5
ind <- c(1,42,80,117,158)
lat_names=rep(NA, 169)
lat_names[ind] <- lat[ind]
lat_names[80] <- 0 #close enough
boxplot(ICV$SST~ICV$Lat, col="green",
#border = rgb(0,0,0,0.3), outpch=20, outcol= rgb(0,0,0,0.3), outbg=rgb(0,0,0,0.3),
ylab="", xlab="", names=lat_names, las=1, cex=1.2)
mtext("SST", side=2, cex=1.5, line=3)
boxplot(ICV$Arag~ICV$Lat, col="green",
ylab="", xlab="", names=lat_names, las=1, cex=1.2)
mtext("Arag.", side=2, cex=1.5, line=3)
boxplot(ICV$pH ~ ICV$Lat, col="green",
ylab="", xlab="", names=lat_names, las=1, cex=1.2)
mtext("pH", side=2, cex=1.5, line=3)
mtext("Latitude", side=1, outer=TRUE, cex=2)
dev.off()
## quartz_off_screen
## 2
##################################
#### Log-transform Arag ####
##################################
dat1$Arag <- log10(dat1$Arag)
dat4.5$Arag <- log10(dat4.5$Arag)
dat8.5$Arag <- log10(dat8.5$Arag)
ICV$Arag <- log10(ICV$Arag)
# Aragonite saturation is a ratio variable; it is limited at zero and
# proportional changes are meaningful. The analysis uses raw values of
# aragonite saturation state. Under raw scaling, the difference between 1 and 2
# is equal to the difference between zero and one, which doesn’t reflect the
# meaning of the variable. Instead, the difference between 1 and 2 should
# be given the same significance as the difference between 0.5 and 1 (doubling vs halving).
# Log-transforming aragonite saturation state will produce this more meaningful
# proportional scaling.
# More generally on the same topic: The need for proportional scaling is the
# reason why precipitation and other ratio variables (e.g. degree-days) are
# typically log-transformed in climate space analysis. In theory, all variables
# should be log-transformed to produce an environmental space where distances
# are comparable for different locations; however, in practice temperature
# doesn’t need to be log-transformed because it doesn’t vary across orders
# of magnitude, and in our case pH is already a log-scaled variable.
#boxplot(dat8.5$pH ~ dat8.5$Year)
#boxplot(dat4.5$pH ~ dat4.5$Year)
# note that these two datasets are exactly the same for years < 2000
#--------------------------------
#### 40-year climate normals ####
#--------------------------------
# for example, 1930 represents from 1/1/1925 to 12/31/1934.
dat_2000 <- dat1 %>% filter(Year>1960 & Year<2010)
dim(dat_2000)
## [1] 1903776 8
# Below is a standardization code that can be deleted becuase the PCA
# does the standardization
# x_SST <- mean(dat_2000$SST)
# s_SST <- sd(dat_2000$SST)
# x_Arag <- mean(dat_2000$Arag)
# s_Arag <- sd(dat_2000$Arag)
# x_pH <- mean(dat_2000$pH)
# s_pH <- sd(dat_2000$pH)
#
# head(dat_2000)
# dat_2000$SST <- (dat_2000$SST - x_SST)/s_SST
# dat_2000$Arag <- (dat_2000$Arag - x_Arag)/s_Arag
# dat_2000$pH <- (dat_2000$pH - x_pH)/s_pH
# dat_2000 %>% summarise_each(sd)
#
dat_1800 <- dat1 %>% filter(Year<1850)
dim(dat_1800)
## [1] 1903776 8
# dat_1800$SST <- (dat_1800$SST - x_SST)/s_SST
# dat_1800$Arag <- (dat_1800$Arag - x_Arag)/s_Arag
# dat_1800$pH <- (dat_1800$pH - x_pH)/s_pH
#
#
dat_2100_8.5 <- dat8.5
# dim(dat_2100_8.5)
# dat_2100_8.5$SST <- (dat_2100_8.5$SST - x_SST)/s_SST
# dat_2100_8.5$Arag <- (dat_2100_8.5$Arag - x_Arag)/s_Arag
# dat_2100_8.5$pH <- (dat_2100_8.5$pH - x_pH)/s_pH
#
dat_2100_4.5 <- dat4.5
# dim(dat_2100_4.5)
# dat_2100_4.5$SST <- (dat_2100_4.5$SST - x_SST)/s_SST
# dat_2100_4.5$Arag <- (dat_2100_4.5$Arag - x_Arag)/s_Arag
# dat_2100_4.5$pH <- (dat_2100_4.5$pH - x_pH)/s_pH
sum(!complete.cases(dat_1800))
## [1] 0
sum(!complete.cases(dat_2000))
## [1] 0
sum(!complete.cases(dat_2100_4.5))
## [1] 0
sum(!complete.cases(dat_2100_8.5))
## [1] 0
sum(!complete.cases(ICV))
## [1] 0
# Calculate climate normals
# Returns the mean climate for each station at each time
norm_1800 <- calculate_normals_annual(dat_1800)
norm_2000 <- calculate_normals_annual(dat_2000)
norm_2100_8.5 <- calculate_normals_annual(dat_2100_8.5)
norm_2100_4.5 <- calculate_normals_annual(dat_2100_4.5)
head(norm_1800)
## No SST Arag pH
## 1 1 -1.2720933 0.2661633 8.129290
## 2 2 -1.1802993 0.2581453 8.136665
## 3 3 -1.0816947 0.2502643 8.141752
## 4 4 -0.9757862 0.2425956 8.144575
## 5 5 -0.8618450 0.2352653 8.145356
## 6 6 -0.7388042 0.2284502 8.144419
dim(norm_1800)
## [1] 39662 4
dim(norm_2000)
## [1] 39662 4
dim(norm_2100_8.5)
## [1] 39662 4
dim(norm_2100_4.5)
## [1] 39662 4
#--------------------------------
### data frame to link station number to Lat Long ####
#--------------------------------
head(dat8.5)
## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 2070 1 0.503480 -0.021016484 7.7819
## 2: 1 20.5 -69.5 2070 2 0.062159 0.003072596 7.8161
## 3: 1 20.5 -69.5 2070 3 -0.320630 -0.015018378 7.8035
## 4: 1 20.5 -69.5 2070 4 -0.688590 -0.024104273 7.8008
## 5: 1 20.5 -69.5 2070 5 -1.676700 -0.042713644 7.7981
## 6: 1 20.5 -69.5 2070 6 0.340380 -0.004470102 7.8048
stations <- unique(dat8.5$No)
stationInfo = data.frame(stations=stations, lat=NA, long=NA)
unik <- which(!duplicated(dat8.5$No))
head(unik)
## [1] 1 49 97 145 193 241
stationInfo$lat <- dat8.5$Lat[unik]
stationInfo$long <- dat8.5$Lon[unik]
names(stationInfo)[1] <- "No"
head(stationInfo)
## No lat long
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5
which(!complete.cases(stationInfo))
## integer(0)
#--------------------------------
### ICV matrix
head(ICV)
## No Lon Lat Year Month SST Arag pH
## 1: 1 20.5 -69.5 1960 1 -0.20420 0.1959550 8.0291
## 2: 1 20.5 -69.5 1960 2 -0.64552 0.2425166 8.0921
## 3: 1 20.5 -69.5 1960 3 -1.02830 0.2176155 8.0698
## 4: 1 20.5 -69.5 1960 4 -1.39630 0.2076613 8.0657
## 5: 1 20.5 -69.5 1960 5 -2.38440 0.1907518 8.0635
## 6: 1 20.5 -69.5 1960 6 -0.36730 0.2253868 8.0699
C <- ICV %>% select(No, SST, Arag, pH)
head(C)
## No SST Arag pH
## 1: 1 -0.20420 0.1959550 8.0291
## 2: 1 -0.64552 0.2425166 8.0921
## 3: 1 -1.02830 0.2176155 8.0698
## 4: 1 -1.39630 0.2076613 8.0657
## 5: 1 -2.38440 0.1907518 8.0635
## 6: 1 -0.36730 0.2253868 8.0699
#--------------------------------
#---------------
### Hemisphere Indexes ###
# N. hemisphere
N_hem <- which(stationInfo$lat>=0)
S_hem <- which(stationInfo$lat<0)
#---------------
#--------------------------------
### 1800 analog to today (Novelty) ####
# What are today's novel climates compared to 1800?
# Novel climates are identified by comparing the
# future climate realization (B) for each gridpoint
# to the past (A) climate realizations (Williams)
# RESTRICT BASELINE TO HEMISPHERE OF OBSERVATION
#--------------------------------
# A is baseline, the base that a station from B is compared to
# If A is past and B is future, then degree of novelty
# If A is future and B is past, then degree of disappearance
# B will be subset to a particular station for analysis
# C is the data frame used to calculate the ICV, will also be subset to a particular station for analysis
A <- norm_1800
# 1800 climate normals
# Serves as the Baseline
head(A)
## No SST Arag pH
## 1 1 -1.2720933 0.2661633 8.129290
## 2 2 -1.1802993 0.2581453 8.136665
## 3 3 -1.0816947 0.2502643 8.141752
## 4 4 -0.9757862 0.2425956 8.144575
## 5 5 -0.8618450 0.2352653 8.145356
## 6 6 -0.7388042 0.2284502 8.144419
dim(A)
## [1] 39662 4
# subset the 2000 data to the same stations as the 1800 data
B <- norm_2000
# 2000 climate normals
head(B)
## No SST Arag pH
## 1 1 -1.4479254 0.1799026 8.035921
## 2 2 -1.3696480 0.1741309 8.045892
## 3 3 -1.2845583 0.1683092 8.053425
## 4 4 -1.1916660 0.1624748 8.058488
## 5 5 -1.0895625 0.1567233 8.061231
## 6 6 -0.9764333 0.1512065 8.062000
dim(B)
## [1] 39662 4
# same columns as A
# sanity check to make sure stations in right order
identical(A$No, B$No) # should be true
## [1] TRUE
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A, B, C, 18906, "_1800_today")
# calc_sigma_D(A, B, C, 18952, "_1800_today")
# calc_sigma_D(A, B, C, 29736, "_1800_today")
# calc_sigma_D(A, B, C, 29789, "_1800_today")
# calc_sigma_D(A, B, C, 8193, "_1800_today")
# calc_sigma_D(A, B, C, 29319, "_1800_today")
#NN.sigma_1800_today <- loop_sigma_D(A, B, C, "_1800_today") # uses whole ocean as a baseline
NN.sigma_1800_today_N <- loop_sigma_D(A[N_hem,], B[N_hem,], C, "_1800_today")
## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413
## [1] 8000 18912 17413
## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_1800_today_S <- loop_sigma_D(A[S_hem,], B[S_hem,], C, "_1800_today")
## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249
## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_1800_today <- rbind(NN.sigma_1800_today_N, NN.sigma_1800_today_S)
head(NN.sigma_1800_today)
## No NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 36 0.033725924 30751 0.23355016
## 2 37 0.015960753 30751 0.16010043
## 3 38 0.010016616 31111 0.12668130
## 4 39 0.001600598 14563 0.05055509
## 5 40 0.018324556 1958 0.17162744
## 6 41 0.004808260 17177 0.08767895
## numPCs_1800_today
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2
hist(NN.sigma_1800_today$numPCs_1800_today)
final_dat <- full_join(stationInfo, NN.sigma_1800_today, by="No")
head(final_dat)
## No lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 1 -69.5 20.5 0.005340623 40066 0.0924152
## 2 2 -68.5 20.5 0.007501424 40009 0.1095737
## 3 3 -67.5 20.5 0.009629160 40067 0.1241975
## 4 4 -66.5 20.5 0.015490226 40009 0.1577081
## 5 5 -65.5 20.5 0.024201364 40066 0.1974682
## 6 6 -64.5 20.5 0.049860481 40066 0.2848815
## numPCs_1800_today
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2
final_dat[which(is.infinite(final_dat$NN.sigma_1800_today)),]
## [1] No lat long
## [4] NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## [7] numPCs_1800_today
## <0 rows> (or 0-length row.names)
final_dat[which(!complete.cases(final_dat)),]
## [1] No lat long
## [4] NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## [7] numPCs_1800_today
## <0 rows> (or 0-length row.names)
#NN.sigma_1800_today[which(is.infinite(NN.sigma_1800_today))] <- NA
identical(nrow(final_dat), nrow(stationInfo))
## [1] TRUE
# should be true, same number as stations as stationInfo
# Visualize
world <- map_data("world2")
Plot_nonInt(final_dat$lat, final_dat$long,
final_dat$NN.sigma_1800_today, world, "sigD novel 2000")
# The NN.sigma_1800_today represents the novelty of today's
# climate compared to 1800
# NN.dist_1800_today represents the Mahalanobis distance from the climate
# in today to it's nearest neighbor in all 1800 data
# No--> NN.station_1800_today represents the station that was queried 2000 (B)
# and it's nearest data in the global reference NN.station_1800_today (1800, A)
# SEE CODE BELOW FOR ADDITION OF NN SITES
#--------------------------------
### Today analog to 1800 Disappearance####
# What are 1800 disappearing climates?
# disappearing climates are identified by comparing
# each past gridpoint (A) to all future climate realizations (B)
#--------------------------------
calc_sigma_D(A, B, C, 1)
## NN.sigma NN.station NN.Mdist num_PCs
## 1 0.005073206 22827 0.09006697 2
calc_sigma_D(B, A, C, 1)
## NN.sigma NN.station NN.Mdist num_PCs
## 1 0.0005525579 38202 0.02969763 2
B <- norm_1800
A <- norm_2000
identical(A$No, B$No) # should be true
## [1] TRUE
identical(sort(A$No), A$No) # should be true
## [1] TRUE
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A, B, C, 18906, "_today_1800")
# calc_sigma_D(A, B, C, 18952, "_today_1800")
# calc_sigma_D(A, B, C, 29736, "_today_1800")
# calc_sigma_D(A, B, C, 29789, "_today_1800")
# calc_sigma_D(A, B, C, 8193, "_today_1800")
# calc_sigma_D(A, B, C, 29319, "_today_1800")
NN.sigma_today_1800_N <- loop_sigma_D(A[N_hem,], B[N_hem,], C, "_today_1800")
## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413
## [1] 8000 18912 17413
## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_today_1800_S <- loop_sigma_D(A[S_hem,], B[S_hem,], C, "_today_1800")
## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249
## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_today_1800 <- rbind(NN.sigma_today_1800_N,
NN.sigma_today_1800_S)
head(NN.sigma_today_1800)
## No NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## 1 36 0.836709022 31931 1.3486467
## 2 37 0.894133247 794 1.4077485
## 3 38 0.661532267 861 1.1633920
## 4 39 0.265489647 39602 0.6854521
## 5 40 0.006257156 39603 0.1000496
## 6 41 0.036266868 39426 0.2423104
## numPCs_today_1800
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2
final_dat2 <- full_join(final_dat, NN.sigma_today_1800)
## Joining, by = "No"
head(final_dat2)
## No lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 1 -69.5 20.5 0.005340623 40066 0.0924152
## 2 2 -68.5 20.5 0.007501424 40009 0.1095737
## 3 3 -67.5 20.5 0.009629160 40067 0.1241975
## 4 4 -66.5 20.5 0.015490226 40009 0.1577081
## 5 5 -65.5 20.5 0.024201364 40066 0.1974682
## 6 6 -64.5 20.5 0.049860481 40066 0.2848815
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.019575390 26505
## 2 2 0.020981196 29981
## 3 2 0.014089519 26082
## 4 2 0.006465863 26082
## 5 2 0.011645698 25749
## 6 2 0.005112661 25973
## NN.Mdist_today_1800 numPCs_today_1800
## 1 0.17743251 2
## 2 0.18374455 2
## 3 0.15036689 2
## 4 0.10170870 2
## 5 0.13663930 2
## 6 0.09041724 2
dim(final_dat2)
## [1] 39662 11
Plot_nonInt(final_dat2$lat, final_dat2$long,
final_dat2$NN.sigma_today_1800, world, "sigD dis 2000")
ggplot(final_dat2) + geom_point(aes(x=NN.sigma_1800_today,
y=NN.sigma_today_1800,
color=lat), alpha=0.5) +
scale_color_gradient2(low="red", high="blue", mid="grey") +
theme_classic() + xlab("Novel climates today from 1800") +
ylab("Disappearing 1800 climates")
# The NN.sigma_1800_disappear represents the novelty of 1800
# compared to today (e.g disappearing climates)
# NN.dist_1800_disappear represents the Mahalanobis distance from the climate
# in 1800 to it's nearest neighbor in all of today's data
# No--> NN.station_1800_disappear represents the the station that was queried 1800 (B)
# and it's nearest data in the global reference NN.station_1800_disappear (2000, A)
stationInfo2 <- stationInfo
names(stationInfo2) <- paste0(names(stationInfo2),"_1800_today")
head(stationInfo2)
## No_1800_today lat_1800_today long_1800_today
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5
dim(stationInfo)
## [1] 39662 3
names(stationInfo2)[1] <-"NN.station_1800_today"
final_dat3 <- left_join(final_dat2, stationInfo2)
## Joining, by = "NN.station_1800_today"
dim(final_dat3)
## [1] 39662 13
head(final_dat3)
## No lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 1 -69.5 20.5 0.005340623 40066 0.0924152
## 2 2 -68.5 20.5 0.007501424 40009 0.1095737
## 3 3 -67.5 20.5 0.009629160 40067 0.1241975
## 4 4 -66.5 20.5 0.015490226 40009 0.1577081
## 5 5 -65.5 20.5 0.024201364 40066 0.1974682
## 6 6 -64.5 20.5 0.049860481 40066 0.2848815
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.019575390 26505
## 2 2 0.020981196 29981
## 3 2 0.014089519 26082
## 4 2 0.006465863 26082
## 5 2 0.011645698 25749
## 6 2 0.005112661 25973
## NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1 0.17743251 2 -69.5 379.5
## 2 0.18374455 2 -69.5 378.5
## 3 0.15036689 2 -68.5 379.5
## 4 0.10170870 2 -69.5 378.5
## 5 0.13663930 2 -69.5 379.5
## 6 0.09041724 2 -69.5 379.5
tail(final_dat3)
## No lat long NN.sigma_1800_today NN.station_1800_today
## 39657 40115 84.5 379.5 0.0033782225 37019
## 39658 40116 85.5 379.5 0.0041649559 37769
## 39659 40117 86.5 379.5 0.0001037297 3229
## 39660 40118 87.5 379.5 0.0047175381 7551
## 39661 40119 88.5 379.5 0.0377153031 7551
## 39662 40120 89.5 379.5 0.0646634655 7551
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 39657 0.07347198 2 1.728258
## 39658 0.08159263 2 1.931682
## 39659 0.01286606 2 2.170640
## 39660 0.08684629 2 2.350300
## 39661 0.24717272 2 2.424026
## 39662 0.32537576 2 2.455741
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800
## 39657 613 2.226041 2
## 39658 613 2.420731 2
## 39659 613 2.648753 2
## 39660 554 2.819972 2
## 39661 3223 2.890216 2
## 39662 3336 2.920432 2
## lat_1800_today long_1800_today
## 39657 86.5 347.5
## 39658 86.5 353.5
## 39659 88.5 59.5
## 39660 88.5 103.5
## 39661 88.5 103.5
## 39662 88.5 103.5
# Should both be 0
sum(!complete.cases(final_dat2))
## [1] 0
sum(!complete.cases(final_dat3))
## [1] 0
head(final_dat3[!complete.cases(final_dat3),])
## [1] No lat long
## [4] NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## [7] numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## [10] NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## [13] long_1800_today
## <0 rows> (or 0-length row.names)
summary(final_dat3, NN.station_today_1800)
## No lat long NN.sigma_1800_today
## Min. : 1 Min. :-78.500 Min. : 20.5 Min. :0.0000002
## 1st Qu.:10100 1st Qu.:-43.500 1st Qu.:138.5 1st Qu.:0.0037303
## Median :20110 Median : -9.500 Median :207.5 Median :0.0778035
## Mean :20095 Mean : -1.731 Mean :208.8 Mean :0.4287141
## 3rd Qu.:30087 3rd Qu.: 34.500 3rd Qu.:295.5 3rd Qu.:0.6417823
## Max. :40120 Max. : 89.500 Max. :379.5 Max. :3.1374782
## NN.station_1800_today NN.Mdist_1800_today numPCs_1800_today
## Min. : 44 Min. :0.000572 Min. :2.000
## 1st Qu.:20665 1st Qu.:0.077212 1st Qu.:2.000
## Median :25736 Median :0.357832 Median :2.000
## Mean :24568 Mean :0.705606 Mean :2.002
## 3rd Qu.:29342 3rd Qu.:1.142396 3rd Qu.:2.000
## Max. :40074 Max. :3.611830 Max. :3.000
## NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## Min. :0.00000 Min. : 39 Min. :0.000314
## 1st Qu.:0.03089 1st Qu.:17767 1st Qu.:0.224749
## Median :0.83348 Median :30424 Median :1.345494
## Mean :0.92699 Mean :24144 Mean :1.296300
## 3rd Qu.:1.60850 3rd Qu.:31698 3rd Qu.:2.112351
## Max. :3.44851 Max. :40119 Max. :3.868075
## numPCs_today_1800 lat_1800_today long_1800_today
## Min. :2.000 Min. :-77.500 Min. : 20.5
## 1st Qu.:2.000 1st Qu.:-43.500 1st Qu.:210.5
## Median :2.000 Median : -2.500 Median :246.5
## Mean :2.002 Mean : -2.377 Mean :244.0
## 3rd Qu.:2.000 3rd Qu.: 25.500 3rd Qu.:284.5
## Max. :3.000 Max. : 89.500 Max. :379.5
final_dat3$latshift_1800_today <- NA
condN <- which(final_dat3$lat>0)
final_dat3$latshift_1800_today[condN] <- final_dat3$lat[condN] -
final_dat3$lat_1800_today[condN]
condS <- which(final_dat3$lat<0)
final_dat3$latshift[condS] <- -1*(final_dat3$lat[condS] -
final_dat3$lat_1800_today[condS])
hist(final_dat3$latshift_1800_today)
#stationInfo3 <- stationInfo
#names(stationInfo3) <- paste0(names(stationInfo3),"_today_1800")
#names(stationInfo3)[1] <- "NN.station_today_1800"
#final_dat4 <- left_join(final_dat3, stationInfo3)
#dim(final_dat4)
#head(final_dat4)
# cond <- which(final_dat4$NN.sigma_1800_today>5 | final_dat4$NN.sigma_1800_disappear>5)
# length(cond)
# final_dat4[cond,]
ggplot(final_dat3) + geom_point(aes(y=lat,
x=lat_1800_today,
color=NN.sigma_1800_today), alpha=0.2) +
scale_color_gradient2(low="red", high="blue", mid="grey") +
theme_classic() + ylab("Latitude in 2000") +
xlab("Latitude of nearest neighbor in 1800") + geom_abline(intercept=0,slope=1)
# This plot makes more sense in the way we typically think about it.
# If we consider a location in 2000, where was the climate it came from in 1800?
write.csv(final_dat3, paste0(results_dir, "1800base_2000_novelty_NShemrestrict.csv"))
final_dat4 <- final_dat3
# ggplot(final_dat4) + geom_point(aes(x=lat_1800_disappear,
# y=lat,
# color=NN.sigma_1800_disappear), alpha=0.2) +
# scale_color_gradient2(low="red", high="blue", mid="grey") +
# theme_classic() + ylab("Latitude in 1800") +
# xlab("Latitude of nearest neighbor in 2000") + geom_abline(intercept=0,slope=1) +
# geom_abline(intercept=0,slope=-1)
# this is weird
#n <- final_dat4 %>% filter(lat_1800_disappear > 0 & lat_1800_disappear < 10)
#Plot_nonInt(n$lat_1800_disappear, n$long_1800_disappear,
# n$NN.sigma_1800_disappear, world, "sigma dis.")
#hist(final_dat4$lat_1800_disappear, breaks=seq(-90,90, by=1))
#hist(final_dat4$lat_1800_today, breaks=seq(-90,90, by=1))
#hist(final_dat4$lat, breaks=seq(-90,90, by=1))
identical(norm_2000$No, norm_2100_8.5$No)
## [1] TRUE
A1 <- norm_2000
# 1970-2000 climate normals
head(A1)
## No SST Arag pH
## 1 1 -1.4479254 0.1799026 8.035921
## 2 2 -1.3696480 0.1741309 8.045892
## 3 3 -1.2845583 0.1683092 8.053425
## 4 4 -1.1916660 0.1624748 8.058488
## 5 5 -1.0895625 0.1567233 8.061231
## 6 6 -0.9764333 0.1512065 8.062000
dim(A1)
## [1] 39662 4
B1 <- norm_2100_8.5
# 2070-2100
head(B1)
## No SST Arag pH
## 1 1 -0.55466154 -0.09210795 7.726612
## 2 2 -0.45603873 -0.09383454 7.740267
## 3 3 -0.34939250 -0.09577437 7.751387
## 4 4 -0.23345842 -0.09793483 7.759896
## 5 5 -0.10654381 -0.10027507 7.765894
## 6 6 0.03343595 -0.10269243 7.769621
dim(B1)
## [1] 39662 4
# same columns as A
# sanity check to make sure stations in right order
identical(A1$No, B1$No) # should be true
## [1] TRUE
identical(A1$No, sort(A1$No)) # should also be true
## [1] TRUE
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A1, B1, C, 18906, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 18952, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 29736, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 29789, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 8193, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C, 29319, "_today_2100_8.5")
# C is the same defined above
NN.sigma_today_2100_8.5_N <- loop_sigma_D(A1[N_hem,], B1[N_hem,], C, "_today_2100_8.5")
## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413
## [1] 8000 18912 17413
## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
# C is the same defined above
NN.sigma_today_2100_8.5_S <- loop_sigma_D(A1[S_hem,], B1[S_hem,], C, "_today_2100_8.5")
## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249
## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_today_2100_8.5 <- rbind(NN.sigma_today_2100_8.5_N,
NN.sigma_today_2100_8.5_S)
which(is.infinite(NN.sigma_today_2100_8.5$NN.sigma_today_2100_8.5))
## integer(0)
which(is.na(NN.sigma_today_2100_8.5$NN.sigma_today_2100_8.5))
## integer(0)
final_dat5 <- full_join(NN.sigma_today_2100_8.5, final_dat4, by="No")
head(final_dat5)
## No NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 36 0.3501681236 29638 0.79989070
## 2 37 0.1043489359 29695 0.41656957
## 3 38 0.0015820951 29638 0.05026184
## 4 39 0.0006816715 28978 0.03298615
## 5 40 0.1188184652 28440 0.44577407
## 6 41 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift
## 1 302.5 26 26
## 2 302.5 27 27
## 3 305.5 -10 -10
## 4 172.5 20 20
## 5 47.5 5 5
## 6 188.5 14 14
dim(final_dat5)
## [1] 39662 19
stationInfo5 <- stationInfo
names(stationInfo5) <- paste0(names(stationInfo),"_today_2100_8.5")
names(stationInfo5)[1] <- "NN.station_today_2100_8.5"
head(stationInfo5)
## NN.station_today_2100_8.5 lat_today_2100_8.5 long_today_2100_8.5
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5
final_dat6 <- left_join(final_dat5, stationInfo5)
## Joining, by = "NN.station_today_2100_8.5"
dim(final_dat6)
## [1] 39662 21
head(final_dat6)
## No NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 36 0.3501681236 29638 0.79989070
## 2 37 0.1043489359 29695 0.41656957
## 3 38 0.0015820951 29638 0.05026184
## 4 39 0.0006816715 28978 0.03298615
## 5 40 0.1188184652 28440 0.44577407
## 6 41 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5
## 1 287.5
## 2 288.5
## 3 287.5
## 4 279.5
## 5 273.5
## 6 271.5
# Visualize
Plot_nonInt(final_dat6$lat, final_dat6$long,
final_dat6$NN.sigma_today_2100_8.5, world, "sigma nov 2100 8.5")
#--------------------------------
### 2100 RCP 8.5 to today, what are today's climates that ####
### will disappear in 2100 RCP 8.5?
#--------------------------------
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D( B1, A1, C, 18906, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 18952, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 29736, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 29789, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 8193, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C, 29319, "_2100_8.5_today")
NN.sigma_2100_8.5_today_N <- loop_sigma_D( B1[N_hem,], A1[N_hem,], C, "_2100_8.5_today")
## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413
## [1] 8000 18912 17413
## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_2100_8.5_today_S <- loop_sigma_D( B1[S_hem,], A1[S_hem,], C, "_2100_8.5_today")
## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249
## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_2100_8.5_today <- rbind(NN.sigma_2100_8.5_today_N, NN.sigma_2100_8.5_today_S)
final_dat7 <- full_join(NN.sigma_2100_8.5_today, final_dat6, by="No")
head(final_dat7)
## No NN.sigma_2100_8.5_today NN.station_2100_8.5_today NN.Mdist_2100_8.5_today
## 1 36 8.292361 40111 10.055106
## 2 37 8.292361 40111 10.340310
## 3 38 8.292361 40111 10.400029
## 4 39 8.292361 40111 10.280321
## 5 40 8.292361 40111 10.091702
## 6 41 8.292361 40111 9.955604
## numPCs_2100_8.5_today NN.sigma_today_2100_8.5 NN.station_today_2100_8.5
## 1 2 0.3501681236 29638
## 2 2 0.1043489359 29695
## 3 2 0.0015820951 29638
## 4 2 0.0006816715 28978
## 5 2 0.1188184652 28440
## 6 2 0.5463419465 28254
## NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5 lat long NN.sigma_1800_today
## 1 0.79989070 2 70.5 20.5 0.033725924
## 2 0.41656957 2 71.5 20.5 0.015960753
## 3 0.05026184 2 72.5 20.5 0.010016616
## 4 0.03298615 2 73.5 20.5 0.001600598
## 5 0.44577407 2 74.5 20.5 0.018324556
## 6 1.03579199 2 75.5 20.5 0.004808260
## NN.station_1800_today NN.Mdist_1800_today numPCs_1800_today
## 1 30751 0.23355016 2
## 2 30751 0.16010043 2
## 3 31111 0.12668130 2
## 4 14563 0.05055509 2
## 5 1958 0.17162744 2
## 6 17177 0.08767895 2
## NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## 1 0.836709022 31931 1.3486467
## 2 0.894133247 794 1.4077485
## 3 0.661532267 861 1.1633920
## 4 0.265489647 39602 0.6854521
## 5 0.006257156 39603 0.1000496
## 6 0.036266868 39426 0.2423104
## numPCs_today_1800 lat_1800_today long_1800_today latshift_1800_today latshift
## 1 2 44.5 302.5 26 26
## 2 2 44.5 302.5 27 27
## 3 2 82.5 305.5 -10 -10
## 4 2 53.5 172.5 20 20
## 5 2 69.5 47.5 5 5
## 6 2 61.5 188.5 14 14
## lat_today_2100_8.5 long_today_2100_8.5
## 1 89.5 287.5
## 2 89.5 288.5
## 3 89.5 287.5
## 4 85.5 279.5
## 5 82.5 273.5
## 6 82.5 271.5
Plot_nonInt(final_dat7$lat, final_dat7$long,
final_dat7$NN.sigma_2100_8.5_today, world, "sigma dis. 2000 in 2100 8.5")
#
# stationInfo7 <- stationInfo
# names(stationInfo7) <- paste0(names(stationInfo),"_2100_8.5_today")
# names(stationInfo7)[1] <- "NN.station_2100_8.5_today"
# head(stationInfo7)
# final_dat8 <- left_join(final_dat7, stationInfo7)
final_dat8 <- final_dat7
#--------------------------------
### 2100 RCP 4.5 against today (novelty) ####
#--------------------------------
identical(norm_2000$No, norm_2100_4.5$No)
## [1] TRUE
# A1 the same as above
B2 <- norm_2100_4.5
# 2070-2100
head(B2)
## No SST Arag pH
## 1 1 -0.8282522 0.04708031 7.879460
## 2 2 -0.7305017 0.04380220 7.891665
## 3 3 -0.6256332 0.04038376 7.901365
## 4 4 -0.5128039 0.03683358 7.908494
## 5 5 -0.3908187 0.03319882 7.913204
## 6 6 -0.2581154 0.02960423 7.915740
dim(B2)
## [1] 39662 4
# same columns as A
# sanity check to make sure stations in right order
identical(A1$No, B2$No) # should be true
## [1] TRUE
identical(A1$No, sort(A1$No)) # should also be true
## [1] TRUE
# C is the same defined above
NN.sigma_today_2100_4.5_N <- loop_sigma_D(A1[N_hem,], B2[N_hem,], C, "_today_2100_4.5")
## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413
## [1] 8000 18912 17413
## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_today_2100_4.5_S <- loop_sigma_D(A1[S_hem,], B2[S_hem,], C, "_today_2100_4.5")
## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249
## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_today_2100_4.5 <- rbind(NN.sigma_today_2100_4.5_N, NN.sigma_today_2100_4.5_S)
which(is.infinite(NN.sigma_today_2100_4.5$NN.sigma_today_2100_4.5))
## integer(0)
which(is.na(NN.sigma_today_2100_4.5$NN.sigma_today_2100_4.5))
## integer(0)
final_dat9 <- full_join(NN.sigma_today_2100_4.5, final_dat8, by="No")
head(final_dat9)
## No NN.sigma_today_2100_4.5 NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5
## 1 36 0.260091777 31342 0.67775158
## 2 37 0.150575943 31338 0.50493740
## 3 38 0.005534165 31456 0.09407847
## 4 39 0.139495959 30347 0.48495881
## 5 40 0.011775140 3437 0.13740012
## 6 41 0.098007656 4227 0.40321232
## numPCs_today_2100_4.5 NN.sigma_2100_8.5_today NN.station_2100_8.5_today
## 1 2 8.292361 40111
## 2 2 8.292361 40111
## 3 2 8.292361 40111
## 4 2 8.292361 40111
## 5 2 8.292361 40111
## 6 2 8.292361 40111
## NN.Mdist_2100_8.5_today numPCs_2100_8.5_today NN.sigma_today_2100_8.5
## 1 10.055106 2 0.3501681236
## 2 10.340310 2 0.1043489359
## 3 10.400029 2 0.0015820951
## 4 10.280321 2 0.0006816715
## 5 10.091702 2 0.1188184652
## 6 9.955604 2 0.5463419465
## NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5 lat
## 1 29638 0.79989070 2 70.5
## 2 29695 0.41656957 2 71.5
## 3 29638 0.05026184 2 72.5
## 4 28978 0.03298615 2 73.5
## 5 28440 0.44577407 2 74.5
## 6 28254 1.03579199 2 75.5
## long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 20.5 0.033725924 30751 0.23355016
## 2 20.5 0.015960753 30751 0.16010043
## 3 20.5 0.010016616 31111 0.12668130
## 4 20.5 0.001600598 14563 0.05055509
## 5 20.5 0.018324556 1958 0.17162744
## 6 20.5 0.004808260 17177 0.08767895
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.836709022 31931
## 2 2 0.894133247 794
## 3 2 0.661532267 861
## 4 2 0.265489647 39602
## 5 2 0.006257156 39603
## 6 2 0.036266868 39426
## NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1 1.3486467 2 44.5 302.5
## 2 1.4077485 2 44.5 302.5
## 3 1.1633920 2 82.5 305.5
## 4 0.6854521 2 53.5 172.5
## 5 0.1000496 2 69.5 47.5
## 6 0.2423104 2 61.5 188.5
## latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 1 26 26 89.5 287.5
## 2 27 27 89.5 288.5
## 3 -10 -10 89.5 287.5
## 4 20 20 85.5 279.5
## 5 5 5 82.5 273.5
## 6 14 14 82.5 271.5
# Visualize
Plot_nonInt(final_dat9$lat, final_dat9$long,
final_dat9$NN.sigma_today_2100_4.5, world, "sig nov 2100 4.5")
stationInfo9 <- stationInfo
names(stationInfo9) <- paste0(names(stationInfo),"_today_2100_4.5")
names(stationInfo9)[1] <- "NN.station_today_2100_4.5"
head(stationInfo9)
## NN.station_today_2100_4.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 1 -69.5 20.5
## 2 2 -68.5 20.5
## 3 3 -67.5 20.5
## 4 4 -66.5 20.5
## 5 5 -65.5 20.5
## 6 6 -64.5 20.5
final_dat10 <- left_join(final_dat9, stationInfo9)
## Joining, by = "NN.station_today_2100_4.5"
dim(final_dat10)
## [1] 39662 31
head(final_dat10)
## No NN.sigma_today_2100_4.5 NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5
## 1 36 0.260091777 31342 0.67775158
## 2 37 0.150575943 31338 0.50493740
## 3 38 0.005534165 31456 0.09407847
## 4 39 0.139495959 30347 0.48495881
## 5 40 0.011775140 3437 0.13740012
## 6 41 0.098007656 4227 0.40321232
## numPCs_today_2100_4.5 NN.sigma_2100_8.5_today NN.station_2100_8.5_today
## 1 2 8.292361 40111
## 2 2 8.292361 40111
## 3 2 8.292361 40111
## 4 2 8.292361 40111
## 5 2 8.292361 40111
## 6 2 8.292361 40111
## NN.Mdist_2100_8.5_today numPCs_2100_8.5_today NN.sigma_today_2100_8.5
## 1 10.055106 2 0.3501681236
## 2 10.340310 2 0.1043489359
## 3 10.400029 2 0.0015820951
## 4 10.280321 2 0.0006816715
## 5 10.091702 2 0.1188184652
## 6 9.955604 2 0.5463419465
## NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5 lat
## 1 29638 0.79989070 2 70.5
## 2 29695 0.41656957 2 71.5
## 3 29638 0.05026184 2 72.5
## 4 28978 0.03298615 2 73.5
## 5 28440 0.44577407 2 74.5
## 6 28254 1.03579199 2 75.5
## long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 20.5 0.033725924 30751 0.23355016
## 2 20.5 0.015960753 30751 0.16010043
## 3 20.5 0.010016616 31111 0.12668130
## 4 20.5 0.001600598 14563 0.05055509
## 5 20.5 0.018324556 1958 0.17162744
## 6 20.5 0.004808260 17177 0.08767895
## numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1 2 0.836709022 31931
## 2 2 0.894133247 794
## 3 2 0.661532267 861
## 4 2 0.265489647 39602
## 5 2 0.006257156 39603
## 6 2 0.036266868 39426
## NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1 1.3486467 2 44.5 302.5
## 2 1.4077485 2 44.5 302.5
## 3 1.1633920 2 82.5 305.5
## 4 0.6854521 2 53.5 172.5
## 5 0.1000496 2 69.5 47.5
## 6 0.2423104 2 61.5 188.5
## latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 1 26 26 89.5 287.5
## 2 27 27 89.5 288.5
## 3 -10 -10 89.5 287.5
## 4 20 20 85.5 279.5
## 5 5 5 82.5 273.5
## 6 14 14 82.5 271.5
## lat_today_2100_4.5 long_today_2100_4.5
## 1 86.5 307.5
## 2 82.5 307.5
## 3 85.5 308.5
## 4 69.5 298.5
## 5 70.5 61.5
## 6 72.5 68.5
#--------------------------------
### today against 2100 RCP 4.5 (disappearing climate) ####
#--------------------------------
NN.sigma_2100_4.5_today_N <- loop_sigma_D(B2[N_hem,], A1[N_hem,], C, "_2100_4.5_today")
## [1] 100 317 17413
## [1] 200 610 17413
## [1] 300 934 17413
## [1] 400 1247 17413
## [1] 500 1683 17413
## [1] 600 2053 17413
## [1] 700 2356 17413
## [1] 800 2660 17413
## [1] 900 2971 17413
## [1] 1000 3213 17413
## [1] 1100 3452 17413
## [1] 1200 3759 17413
## [1] 1300 3996 17413
## [1] 1400 4236 17413
## [1] 1500 4540 17413
## [1] 1600 4847 17413
## [1] 1700 5220 17413
## [1] 1800 5530 17413
## [1] 1900 5830 17413
## [1] 2000 6129 17413
## [1] 2100 6427 17413
## [1] 2200 6726 17413
## [1] 2300 7099 17413
## [1] 2400 7694 17413
## [1] 2500 8267 17413
## [1] 2600 8629 17413
## [1] 2700 8887 17413
## [1] 2800 9069 17413
## [1] 2900 9286 17413
## [1] 3000 9454 17413
## [1] 3100 9622 17413
## [1] 3200 9790 17413
## [1] 3300 9954 17413
## [1] 3400 10119 17413
## [1] 3500 10287 17413
## [1] 3600 10458 17413
## [1] 3700 10592 17413
## [1] 3800 10764 17413
## [1] 3900 10945 17413
## [1] 4000 11095 17413
## [1] 4100 11309 17413
## [1] 4200 11469 17413
## [1] 4300 11702 17413
## [1] 4400 11871 17413
## [1] 4500 12109 17413
## [1] 4600 12278 17413
## [1] 4700 12450 17413
## [1] 4800 12692 17413
## [1] 4900 12863 17413
## [1] 5000 13035 17413
## [1] 5100 13280 17413
## [1] 5200 13453 17413
## [1] 5300 13627 17413
## [1] 5400 13878 17413
## [1] 5500 14056 17413
## [1] 5600 14235 17413
## [1] 5700 14414 17413
## [1] 5800 14670 17413
## [1] 5900 14848 17413
## [1] 6000 15026 17413
## [1] 6100 15204 17413
## [1] 6200 15384 17413
## [1] 6300 15641 17413
## [1] 6400 15819 17413
## [1] 6500 15998 17413
## [1] 6600 16176 17413
## [1] 6700 16356 17413
## [1] 6800 16534 17413
## [1] 6900 16791 17413
## [1] 7000 16970 17413
## [1] 7100 17149 17413
## [1] 7200 17328 17413
## [1] 7300 17506 17413
## [1] 7400 17684 17413
## [1] 7500 17862 17413
## [1] 7600 18119 17413
## [1] 7700 18297 17413
## [1] 7800 18476 17413
## [1] 7900 18655 17413
## [1] 8000 18912 17413
## [1] 8100 19090 17413
## [1] 8200 19269 17413
## [1] 8300 19526 17413
## [1] 8400 19704 17413
## [1] 8500 19881 17413
## [1] 8600 20058 17413
## [1] 8700 20312 17413
## [1] 8800 20488 17413
## [1] 8900 20664 17413
## [1] 9000 20916 17413
## [1] 9100 21091 17413
## [1] 9200 21266 17413
## [1] 9300 21441 17413
## [1] 9400 21691 17413
## [1] 9500 21866 17413
## [1] 9600 22041 17413
## [1] 9700 22218 17413
## [1] 9800 22468 17413
## [1] 9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_2100_4.5_today_S <- loop_sigma_D(B2[S_hem,], A1[S_hem,], C, "_2100_4.5_today")
## [1] 100 140 22249
## [1] 200 298 22249
## [1] 300 459 22249
## [1] 400 622 22249
## [1] 500 765 22249
## [1] 600 909 22249
## [1] 700 1052 22249
## [1] 800 1195 22249
## [1] 900 1345 22249
## [1] 1000 1468 22249
## [1] 1100 1619 22249
## [1] 1200 1763 22249
## [1] 1300 1923 22249
## [1] 1400 2093 22249
## [1] 1500 2228 22249
## [1] 1600 2401 22249
## [1] 1700 2539 22249
## [1] 1800 2718 22249
## [1] 1900 2858 22249
## [1] 2000 3040 22249
## [1] 2100 3182 22249
## [1] 2200 3372 22249
## [1] 2300 3518 22249
## [1] 2400 3710 22249
## [1] 2500 3856 22249
## [1] 2600 4047 22249
## [1] 2700 4190 22249
## [1] 2800 4368 22249
## [1] 2900 4506 22249
## [1] 3000 4678 22249
## [1] 3100 4810 22249
## [1] 3200 4967 22249
## [1] 3300 5092 22249
## [1] 3400 5240 22249
## [1] 3500 5367 22249
## [1] 3600 5531 22249
## [1] 3700 5666 22249
## [1] 3800 5834 22249
## [1] 3900 5969 22249
## [1] 4000 6141 22249
## [1] 4100 6278 22249
## [1] 4200 6451 22249
## [1] 4300 6586 22249
## [1] 4400 6753 22249
## [1] 4500 6882 22249
## [1] 4600 7031 22249
## [1] 4700 7149 22249
## [1] 4800 7275 22249
## [1] 4900 7388 22249
## [1] 5000 7514 22249
## [1] 5100 7640 22249
## [1] 5200 7753 22249
## [1] 5300 7882 22249
## [1] 5400 8010 22249
## [1] 5500 8139 22249
## [1] 5600 8256 22249
## [1] 5700 8404 22249
## [1] 5800 8538 22249
## [1] 5900 8680 22249
## [1] 6000 8864 22249
## [1] 6100 9073 22249
## [1] 6200 9259 22249
## [1] 6300 9501 22249
## [1] 6400 9758 22249
## [1] 6500 10019 22249
## [1] 6600 10336 22249
## [1] 6700 10614 22249
## [1] 6800 10899 22249
## [1] 6900 11128 22249
## [1] 7000 11363 22249
## [1] 7100 11532 22249
## [1] 7200 11769 22249
## [1] 7300 11938 22249
## [1] 7400 12178 22249
## [1] 7500 12350 22249
## [1] 7600 12523 22249
## [1] 7700 12772 22249
## [1] 7800 12948 22249
## [1] 7900 13203 22249
## [1] 8000 13382 22249
## [1] 8100 13561 22249
## [1] 8200 13824 22249
## [1] 8300 14012 22249
## [1] 8400 14279 22249
## [1] 8500 14462 22249
## [1] 8600 14649 22249
## [1] 8700 14915 22249
## [1] 8800 15102 22249
## [1] 8900 15289 22249
## [1] 9000 15556 22249
## [1] 9100 15740 22249
## [1] 9200 15926 22249
## [1] 9300 16112 22249
## [1] 9400 16386 22249
## [1] 9500 16573 22249
## [1] 9600 16760 22249
## [1] 9700 16947 22249
## [1] 9800 17224 22249
## [1] 9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_2100_4.5_today <- rbind(NN.sigma_2100_4.5_today_N, NN.sigma_2100_4.5_today_S)
which(is.infinite(NN.sigma_2100_4.5_today$NN.sigma_2100_4.5_today))
## integer(0)
which(is.na(NN.sigma_2100_4.5_today$NN.sigma_2100_4.5_today))
## integer(0)
final_dat11 <- full_join(NN.sigma_2100_4.5_today, final_dat10, by="No")
head(final_dat11)
## No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36 3.551980 40055 3.967162
## 2 37 3.505984 40055 3.923107
## 3 38 3.315848 40055 3.741139
## 4 39 2.951652 40055 3.393209
## 5 40 2.553620 40111 3.013691
## 6 41 2.223800 40111 2.699426
## numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1 2 0.260091777 31342
## 2 2 0.150575943 31338
## 3 2 0.005534165 31456
## 4 2 0.139495959 30347
## 5 2 0.011775140 3437
## 6 2 0.098007656 4227
## NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1 0.67775158 2 8.292361
## 2 0.50493740 2 8.292361
## 3 0.09407847 2 8.292361
## 4 0.48495881 2 8.292361
## 5 0.13740012 2 8.292361
## 6 0.40321232 2 8.292361
## NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1 40111 10.055106 2
## 2 40111 10.340310 2
## 3 40111 10.400029 2
## 4 40111 10.280321 2
## 5 40111 10.091702 2
## 6 40111 9.955604 2
## NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 0.3501681236 29638 0.79989070
## 2 0.1043489359 29695 0.41656957
## 3 0.0015820951 29638 0.05026184
## 4 0.0006816715 28978 0.03298615
## 5 0.1188184652 28440 0.44577407
## 6 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 287.5 86.5 307.5
## 2 288.5 82.5 307.5
## 3 287.5 85.5 308.5
## 4 279.5 69.5 298.5
## 5 273.5 70.5 61.5
## 6 271.5 72.5 68.5
# Visualize
Plot_nonInt(final_dat11$lat, final_dat11$long,
final_dat11$NN.sigma_2100_4.5_today, world, "sigma dis. 2000 in 2100 4.5")
# stationInfo11 <- stationInfo
# names(stationInfo11) <- paste0(names(stationInfo),"_2100_4.5_today")
# names(stationInfo11)[1] <- "NN.station_2100_4.5_today"
# head(stationInfo11)
# final_dat12 <- left_join(final_dat11, stationInfo11)
final_dat12 <- final_dat11
dim(final_dat12)
## [1] 39662 35
head(final_dat12)
## No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36 3.551980 40055 3.967162
## 2 37 3.505984 40055 3.923107
## 3 38 3.315848 40055 3.741139
## 4 39 2.951652 40055 3.393209
## 5 40 2.553620 40111 3.013691
## 6 41 2.223800 40111 2.699426
## numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1 2 0.260091777 31342
## 2 2 0.150575943 31338
## 3 2 0.005534165 31456
## 4 2 0.139495959 30347
## 5 2 0.011775140 3437
## 6 2 0.098007656 4227
## NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1 0.67775158 2 8.292361
## 2 0.50493740 2 8.292361
## 3 0.09407847 2 8.292361
## 4 0.48495881 2 8.292361
## 5 0.13740012 2 8.292361
## 6 0.40321232 2 8.292361
## NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1 40111 10.055106 2
## 2 40111 10.340310 2
## 3 40111 10.400029 2
## 4 40111 10.280321 2
## 5 40111 10.091702 2
## 6 40111 9.955604 2
## NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 0.3501681236 29638 0.79989070
## 2 0.1043489359 29695 0.41656957
## 3 0.0015820951 29638 0.05026184
## 4 0.0006816715 28978 0.03298615
## 5 0.1188184652 28440 0.44577407
## 6 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 287.5 86.5 307.5
## 2 288.5 82.5 307.5
## 3 287.5 85.5 308.5
## 4 279.5 69.5 298.5
## 5 273.5 70.5 61.5
## 6 271.5 72.5 68.5
sum(!complete.cases(final_dat12))
## [1] 22249
### Final dataframe ####
dat_Nov <- final_dat12
### Write to file ####
write.csv(final_dat12, paste0(results_dir, "SigmaD_Hem.csv"), row.names = FALSE)
head(dat_Nov)
## No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36 3.551980 40055 3.967162
## 2 37 3.505984 40055 3.923107
## 3 38 3.315848 40055 3.741139
## 4 39 2.951652 40055 3.393209
## 5 40 2.553620 40111 3.013691
## 6 41 2.223800 40111 2.699426
## numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1 2 0.260091777 31342
## 2 2 0.150575943 31338
## 3 2 0.005534165 31456
## 4 2 0.139495959 30347
## 5 2 0.011775140 3437
## 6 2 0.098007656 4227
## NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1 0.67775158 2 8.292361
## 2 0.50493740 2 8.292361
## 3 0.09407847 2 8.292361
## 4 0.48495881 2 8.292361
## 5 0.13740012 2 8.292361
## 6 0.40321232 2 8.292361
## NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1 40111 10.055106 2
## 2 40111 10.340310 2
## 3 40111 10.400029 2
## 4 40111 10.280321 2
## 5 40111 10.091702 2
## 6 40111 9.955604 2
## NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 0.3501681236 29638 0.79989070
## 2 0.1043489359 29695 0.41656957
## 3 0.0015820951 29638 0.05026184
## 4 0.0006816715 28978 0.03298615
## 5 0.1188184652 28440 0.44577407
## 6 0.5463419465 28254 1.03579199
## numPCs_today_2100_8.5 lat long NN.sigma_1800_today NN.station_1800_today
## 1 2 70.5 20.5 0.033725924 30751
## 2 2 71.5 20.5 0.015960753 30751
## 3 2 72.5 20.5 0.010016616 31111
## 4 2 73.5 20.5 0.001600598 14563
## 5 2 74.5 20.5 0.018324556 1958
## 6 2 75.5 20.5 0.004808260 17177
## NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1 0.23355016 2 0.836709022
## 2 0.16010043 2 0.894133247
## 3 0.12668130 2 0.661532267
## 4 0.05055509 2 0.265489647
## 5 0.17162744 2 0.006257156
## 6 0.08767895 2 0.036266868
## NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1 31931 1.3486467 2 44.5
## 2 794 1.4077485 2 44.5
## 3 861 1.1633920 2 82.5
## 4 39602 0.6854521 2 53.5
## 5 39603 0.1000496 2 69.5
## 6 39426 0.2423104 2 61.5
## long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1 302.5 26 26 89.5
## 2 302.5 27 27 89.5
## 3 305.5 -10 -10 89.5
## 4 172.5 20 20 85.5
## 5 47.5 5 5 82.5
## 6 188.5 14 14 82.5
## long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1 287.5 86.5 307.5
## 2 288.5 82.5 307.5
## 3 287.5 85.5 308.5
## 4 279.5 69.5 298.5
## 5 273.5 70.5 61.5
## 6 271.5 72.5 68.5
for_liqing <- data.frame(No=dat_Nov$No, lat=dat_Nov$lat,long=dat_Nov$long,
A=dat_Nov$NN.sigma_today_1800,
B=dat_Nov$NN.sigma_1800_today,
C=dat_Nov$NN.sigma_2100_4.5_today,
D=dat_Nov$NN.sigma_today_2100_4.5,
E=dat_Nov$NN.sigma_2100_8.5_today,
F=dat_Nov$NN.sigma_today_2100_8.5)
write.csv(for_liqing, paste0(results_dir,"SigmaD_Hem_plot1.csv"), row.names=FALSE)
for_liqing_Md <- data.frame(No=dat_Nov$No, lat=dat_Nov$lat,long=dat_Nov$long,A=dat_Nov$NN.Mdist_today_1800, B=dat_Nov$NN.Mdist_1800_today, C=dat_Nov$NN.Mdist_2100_4.5_today, D=dat_Nov$NN.Mdist_today_2100_4.5, E=dat_Nov$NN.Mdist_2100_8.5_today, F=dat_Nov$NN.Mdist_today_2100_8.5)
write.csv(for_liqing_Md, paste0(results_dir,"SigmaD_Hem_plot2.csv"), row.names=FALSE)
#plot(final_dat12$lat_today_2100_8.5, final_dat12$lat) #If A is today and B is future,
# where station No's climate in the future will come from today (or the closest similar climate today).
# Latitude of the station today (x-axis) where the latitude of the station in the future (y-axis) will come from
#plot(final_dat12$lat, final_dat12$lat_2100_8.5_today) # If A is future and B is today, this
# represents where station "No" will be found in the future (or the closest similar climate).
# Latitude of the station today (x-axis) and where it will be in the future (y-axis)
### Calculate percentages ####
N <- nrow(dat_Nov)
calc_perc <- function(x, descrip){
print(paste("Percent of cells with moderate dissimilarity for", descrip,":",
round(sum(x>2 & x<4)/N,3)*100, sep=" "))
print(paste("Percent of cells with extreme dissimilarity for", descrip,":",
round(sum(x>4)/N, 3)*100, sep=" "))
}
calc_perc(dat_Nov$NN.sigma_today_1800, "1800 disappearing in 2000")
## [1] "Percent of cells with moderate dissimilarity for 1800 disappearing in 2000 : 12.4"
## [1] "Percent of cells with extreme dissimilarity for 1800 disappearing in 2000 : 0"
calc_perc(dat_Nov$NN.sigma_1800_today, "Novel climates in 2000 since 1800")
## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2000 since 1800 : 3.7"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2000 since 1800 : 0"
calc_perc(dat_Nov$NN.sigma_2100_4.5_today, "2000 disappearing in 2100 4.5")
## [1] "Percent of cells with moderate dissimilarity for 2000 disappearing in 2100 4.5 : 34.2"
## [1] "Percent of cells with extreme dissimilarity for 2000 disappearing in 2100 4.5 : 35.6"
calc_perc(dat_Nov$NN.sigma_today_2100_4.5, "Novel climates in 2100 4.5 since 2000")
## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2100 4.5 since 2000 : 42.1"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2100 4.5 since 2000 : 10.3"
calc_perc(dat_Nov$NN.sigma_2100_8.5_today, "2000 disappearing in 2100 8.5")
## [1] "Percent of cells with moderate dissimilarity for 2000 disappearing in 2100 8.5 : 3.5"
## [1] "Percent of cells with extreme dissimilarity for 2000 disappearing in 2100 8.5 : 95"
calc_perc(dat_Nov$NN.sigma_today_2100_8.5, "Novel climates in 2100 8.5 since 2000")
## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2100 8.5 since 2000 : 6.7"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2100 8.5 since 2000 : 81.9"
# Visualize relationship between sigma_D and M_D
plot(dat_Nov$NN.sigma_today_1800, dat_Nov$NN.Mdist_today_1800)
plot(dat_Nov$NN.sigma_today_2100_8.5, dat_Nov$NN.Mdist_today_2100_8.5, col=dat_Nov$numPCs_today_2100_8.5)
legend(0, 15, 1:4, 1:4)
abline(h=12)
hist(dat_Nov$numPCs_today_2100_8.5)